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in Incompressible Boundary Layers 


P. Balakumar 

Aerospace Engineering Department, 

Old Dominion University, Norfolk, Virginia 23529-0247 

Abstract 

We have developed a code where the nonlinear terms are treated implicitly. The equations 
are discretized using the two-point fourth order compact scheme in the y-direction and the 
backward Euler method in the x-direction. We investigated the transition process in a Blasius 
boundary layer due to fundamental type breakdown. With 8 modes in the u> and 3 planes, we 
could compute the evolution of disturbances up to Re x =910, which is well into the strongly 
nonlinear region. The transition onset point is located around Re x =850. The comparison with 
the measurements and with the DNS computations are very good up to Re x =880. 

1. Introduction 

Breakdown from laminar to turbulent flow in zero and mild pressure gradient boundary layers 
is caused by Tollmien-Schlichting (viscous) instability. Though there exist several mechanisms 
and routes to go from a laminar to a turbulent state, all of them in general follow these 
fundamental processes: 

Receptivity 

Linear instability 

Nonlinear stability and saturation 

Secondary instability 

Breakdown to turbulence 
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In the receptivity process, the unsteadiness in the environment and the inhomogeneities in the 
geometry generate instability waves inside the flow. In quiet environments, the initial amplitudes 
of these instability waves are small compared to any characteristic velocity and length scales in 
the flow. Goldstein (1983 a) theoretically explained using asymptotic methods how the Tollmien- 
Schlichting waves (T-S waves ) are generated near a leading edge of a flat plate by the long 
wavelength acoustic disturbances and in a companion paper (1985) described the scattering of 
T-S waves from the acoustic disturbances by the streamwise variations in surface geometries. 
In the second stage, the amplitudes of these instability waves grow exponentially downstream 
and this process is governed by the linearized Navier-Stokes equations. Further downstream, 
the amplitudes of the disturbances become large and the nonlinear effects inhibit the exponential 
growth and the amplitudes of the disturbances eventually saturate or attain singular values. In the 
next stage, these finite amplitude saturated disturbances become unstable to two- and/or three- 
dimensional disturbances. This is called secondary instability and beyond this stage the spectrum 
broadens, due to complex interactions and further instabilities, and the flow becomes turbulent 
in a short distance downstream. In this paper, we investigated the later stages of transition in 
Blasius boundary layers using Parabolized Stability Equations ( PSE ) approach. 

The Parabolized Stability Equations (PSE) approach currently can predict the first three 
stages of the transition process, linear instability, nonlinear stability and saturation and secondary 
instability, accurately and very efficiently, Herbert (1997). After the skin friction rise the PSE 
computations cease to converge. There may be two sources for the program not to converge in 
this region. One is that the PSE approximation itself may not be valid in the highly nonlinear 
region. The other may be that the iteration on the nonlinear terms, at present they are treated 
explicitly and iterated till they converge, may not be converging. In this work, we will treat 
the nonlinear terms implicitly and will investigate how far downstream we can continue the 
PSE computations. If we can compute up to the skin friction maximum in the transition region, 
this will help us in developing new transition modeling and in developing improved transition 
prediction methods. We also observed that in some PSE computations, e.g., Gortler instability, 
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crossflow instability, at the later stages, the meanflow distortion term converges very slowly. 
Hence, in these cases the implicit formulation will improve the convergence. 


2. Formulation 


In the parabolized stability equations (PSE) approach, one attempts to construct an approxi- 
mate solution of the full Navier-Stokes equations. The idea was first introduced by Herbert( 1991 ) 
and applied to linear and non-linear Blasius boundary-layer flow by Bertolotti (1992). Now it 
has been developed and has been applied to two and three-dimensional incompressible and com- 
pressible boundary-layer flows (Chang et.al. 1994, Malik et.al 1994). Herbert (1997) in a recent 
review described the development and the application of PSE to different problems and here we 
give the governing equations and the procedure that we use to solve the equations. 

We investigate the transition process in a Blasius boundary layer over a flat plate. Let us 
denote the Cartesian coordinate system by x* (i= 1,2,3), the velocity components by u, and the 
pressure by p. We decompose the total flow quantities as the sum of the mean flow and the 
disturbance quantities. 


«,(x, t) = Uoi(x) + u,(x,t ), 


Pi(x,t) = P 0 (x) -f Pi(x,t), 


( 1 ) 


where, Uoi are the mean velocity components which are the solution of the Blasius equation 
and P is the mean pressure. Substituting these expressions into the Navier-Stokes equations we 
obtain the nonlinear equations for the disturbances in the form 


du x 

dx x 


= o. 


dui dUoi 

~m +Ul ~d. 


+ 


du, du t 
° J d.Vj + Uj dxj 



i d 2 u, 
Re dxjdxj 


( 2 ) 


where the Reynolds number 


Re = 


UooL 

1 

V 


(3) 
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, Uoc is the freestream velocity, and L is a reference length. 

To derive the PSE equations, we write the dependent variables «i, i/>, 113 , p as the sum of 
normal modes 


MB 


MO 


- 2! X! qmn 

m=-MB n=-MO 


> \l / »mn(Xl )dx \+i7n3£2-inujt 


( 4 ) 


Here u is the frequency of the disturbance, j3 is the spanwise wave number, and Q mn (xi) is the 
wave number in the axial direction. q mn (xi ,X 3 ) is the amplitude function for the mode (m,n) 
, MO is the total number of modes kept in the frequency domain in one quadrant and MB is 
the number of spanwise modes. Substituting this expression into the disturbance equations and 
making the PSE approximation, we obtain PSE equations for each mode (m,n). 


du3 mn . ~ _ , ■ j _ 

-f* lOtjjinUlmn ”1“ ITUfjUo mn i 

OX 3 


dll Imn 
dx i 


= 0. 


( 5 ) 
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T Re— U 3 mn T l&mnPmn 
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dx 1 
dui 




diii 
dx 3 
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( 8 ) 


At present, the nonlinear terms S mn are evaluated from the values obtained at the previous 
iteration step and using the FFT technique. In the implicit formulation, these terms are function 
of Q mn and the system of equations becomes a set of coupled nonlinear equations. This nonlinear 
coupled equations will be solved using Newton - Raphson iteration technique Balakumar (1997). 


3. Solution procedure 

The eqs. (5-8) are solved using the two-point fourth order compact scheme, Malik et. al 
(1982), derived by means of the Euler-Maclaurin formula 


0* 


0*" 1 


h k (dil> k d 
2 \dY + dY ) 


+0(hl ), 


h\ ( d 2 p k d 2 ib k 1 \ 
12 ( dY 1 + dY 2 ) 


(9) 


where, 0 k = 0(Vfc) and h k - Y k - Y k _\. 

To apply the above scheme to this problem, we have to rewrite eqs. (5-8) as a system of 
first order equations. The eqs. (5-8) can easily be written as 


0 * 


= F n .m(4\ Y) n = -MB to MB, m = -MO to MO, 


(10) 
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where, p = ^ A/B,— MO' ^’A/fl,A/o} a nd 4'n.in = { u n,mi w n.mi w n,m’’ v n,m,Pn,m } • 

4' n ,m and Fn,m are (6x1) complex column vectors and 4' is a (6 x ( 2MB + \)*(2MO + 1)) 
complex column vector. From eq. (10), we get 


" 

v *- m dY + 


f dF n ,, 

1 dp 


F, 


(ID 


where, F = {^Tv/B -a/O’ •••’ ^A/B MOj- Substituting eqs. (10) and (1 1) into (9), we obtain 


*lm - rf:™ = + f ‘™‘) 

« ’ 


' t | w/t , i' dF n< 

*- n m i 


rt _ piik — 1 
1 n,m 


f>p \ k — i 
r n.m \ pk~ 1 


( 12 ) 


12 ^ n ’ m \ dp J n m V dp 

for n = — M B to MB, m = —MO to MO and A: = 2, )V + l. 


where, (V + 1 is the total number of grid points in the normal direction. This is a system of 
non-linear algebraic equations which is solved using the Newton-Raphson iterative method. We 
write pn,m = Pon.m + AV>»,m and substitute into the eq. (12) and linearize in Ap, to obtain 


Pon.m + A Pn,m ~ Pon,m ~ &Pn,m 
^k — 1 


h ± F * 

r> 1 on.m 1 ^ on,m 
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+ FL~L + [^- ) AP k +' 9Fc 
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hi 
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dip 
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on ’ m V dp J 
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dp J \ t 


f) 2 p 

k , / t/ -T on,m 


dp' 1 


k - 1 


- 
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d/; 


fc-i 

on.m \ zrt-1 


(13) 


dp 


Ap 


k - 1 


d^Fon 

dp 2 


k - 1 


F k ~ l Ap k ~ 1 1 

for « = 0, M B. m = 0, MO, and A - = 2, iV + 1. 
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Collecting the terms we obtain a linear system of equations for Ac* and A c* 1 which can 
be written as 


A k &4’ k ~ l + B k Ap k = D k for k = 2, N + 1, 


(14) 


where .4*, B k are complex matrices of size ((MB + l)x(MO+l) x 6, (2MB + l)x(2MO+l) x 6) and 
D k is a complex vector of size (MB + l)x(MO+l) x 6. Using the property A 0„, m = At _ m 
A i'n,m = Ac'-;,. m (where * stands for complex conjugate), we can rewrite this equation in real 
and imaginary parts of Ac„, m (n = 0 m = 0, ...,MO). 


A k 


A t/v 
Ac, 




A: = 2, ( A r + 1), 


(15) 


where /1/t, B k are real matrices of size (MB + l)x(MO+l) x 12. The boundary conditions can 
be written as. 


fit 



= fii, 


(16) 


and 


'4/V+2 


Ai/v 

A0, 


,v+i 

= fi.V + 2. 


(17) 


where, and .4^+2 are ((MB + l)x(MO+l) x 6, (MB + l)x(MO+l) x 12) real matrices and 
Du fijV+2 are ((MB + l)x(MO+l) x 6) column vectors. Combining eqs. (15-17) we can form 
a block tridiagonal system or can write it as a banded system, both of which can be solved 
very efficiently. 
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4. Computational efficiency 


In order to discuss the efficiency and the robustness of the implicit solution method in 
comparison with the explicit method, this section will show results from an application of both 
methods to the Swept Hiemenz flow problem. Relevant parameters and the mean flow profiles 
are given in Malik, Fei and Chang (1996). For a Reynolds number in the spanwise direction 
of Re o = 500. a steady disturbance with a spanwise wave number of d\ = 0.4 and an initial 
amplitude of A=0 . \% is introduced at the streamwise position of Re o = 150. 

Figure 1 shows the growth rates a\ based on the disturbance component in the spanwise 
direction that are obtained from the explicit and the implicit method. Additionally, the linear 
growth rates are plotted. Besides the perfect agreement of the results obtained from the different 
methods, it is noted that the nonlinear effects modify the growth rate starting at a Reynolds 
number of Re=400. It is further seen that the explicit method ceases to converge at Re=570 for 
an SUR-parameter (Successive Under Relaxation) of 1 .0., whereas the implicit method continues 
to converge well beyond the plotted region. It is noted that a Reynolds number of Re=570 
corresponds to the region where the stationary disturbances saturate, and hence, the nonlinear 
effects are indeed strong at the point where the explicit method stops converging for an SUR- 
parameter of 1.0. 

The number of iterations on the nonlinear terms necessary to obtain converged solution, and 
the convergence history at Re=570 are shown in figures 2 and 3, respectively. There, the residue 
is defined as the maximal difference in the shape functions of the primary disturbance obtained in 
two consecutive iterations. The extremely rapid convergence of the implicit method is visualized 
in the figure 3 and shows the superiority of the implicit method in the developed nonlinear region. 

Two observations can be made from these plots. First, in the region of linear disturbance 
growth (Re<400), the explicit method needs as many iterations on the nonlinear terms as the 
implicit method. Since the explicit method iterates just on the primary disturbance, whereas the 
implicit method computes all disturbances simultaneously, the explicit method is expected to be 
much more efficient in that region. With an increasing effect of the nonlinear terms for Re>400, 
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however, this similar performance of the two methods is altered significantly. The number of 
iterations on the nonlinear terms using the explicit method increases exponentially, whereas the 
performance of the implicit code basically remains unchanged. 

It is concluded that the implicit method provides a powerful tool to conduct nonlinear PSE 
computations in the regions of highly nonlinear disturbance growth. The most effective technique 
in applying the developed algorithms is seen in combining the explicit and implicit methods by 
starting the computations with the explicit method, and once a critical number of nonlinear 
iterations is reached, by continuing with the implicit method. 


5. Transition due to Klebanoff type breakdown 

Next we investigated the transition process in the Blasius boundary layer due to fundamental 
breakdown process (K-type). We simulated one of the experiments performed by Saric et al. 
(1981) and the parameters are given in Table 1. 

F = .76* 10 -4 
13 = .142 

/?e 0 = 710 (18) 

( u 0,l )max = 00:338 

= -Q 0024 


Here, the variables are nondimensionalised by 


velocity : Free stream velocity Uoo 


length : 



Here, F is the non-dimensional frequency frequency defined by 


(19) 


F = 


'Ixvf 
U- 1 

^ oo 


( 20 ) 
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and Reo is the Reynolds number at the initial location x=xo- The Reynolds numbers Re x , RE* 
and the skin friction coefficient Cf are defined by 


Re x = 
RE; 

Cf = 


too* 

v 

L'ooX 

V 


pUl 


( 21 ) 


Figure 4 shows the variation of the maximum rms fluctuations for the streamwise velocity 
with the Reynolds number. The results are compared with the experimental results of Saric et. 
al. (1984) and with the DNS results of Liu et. al. (1995). We used 8 and 12 modes, MB=MO=8, 
12 in these computations. With the SUR parameter ( Successive under relaxation) equals to 1.0, 
the explicit code ceases to converge beyond Re x =860, and with SUR=.25, we can extend the 
computations up to Re x =880 and with the implicit code up to Re x =905 with 8 modes. With the 
12 modes, the explicit code stops converging at Re x =880 and we did not have enough memory 
in the computer to continue further using the implicit code. It is seen that the computed results 
are in good agreement with the experiment and with the DNS results. 

Figure 5 shows the distribution of u^ in the y-direction at different axial locations Re x =855, 
865, 875, 885 and 890. Re x =855 is approximately the transition onset point. We observe that the 
location where u rms peaks increases with the downstream direction. This is in good qualitative 
agreement with the Klebanoff et. al. (1961) measurements. 

Figure 6 shows the variation of the skin friction coefficient Cf with the Reynolds number. For 
comparison, we also plotted the skin friction coefficient for the Blasius boundary layer. We note 
that the transition onset point is located at Re x =855 and beyond this point skin friction increases 
sharply. These results show that we can compute the transition process in the fully nonlinear 
region using PSE approach. In figure 7, we plotted the modified mean velocity profile and the 
unmodified Blasius boundary layer profile at Re x =890. It is seen that the velocity increases in 
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the region close to the wall and decreases in the outer part of the boundary layer. Hence the 
boundary layer profiles gradually become fuller in the fully nonlinear region. 

In figure 8, we plotted the variation of unsteady u fluctuations at the maximum rms locations 
with time at different Reynolds numbers. The results show that up to Re x =875, the oscillations are 
approximately sinusoidal and beyond this stage the disturbances have smaller positive amplitude 
and larger negative value. At Re x =890, the maximum positive value is about 0.1 and the 
minimum value is -0.2. Hence it shows that the evolution of the disturbances approach 1 
“spike” stage according to Klebanoff et. al. (1961). 

In figures 9(a-i) we plotted the evolution of different modes in the streamwise direction. 
Each figure shows the results for at a fixed mw for different nfi where n=0 to 8. Figure 9(a) 
shows the results for MO=0 and figure 9(i) depicts the results for MO=9. First observation is that 
the amplitude of the disturbances gradually decreases with increasing mw. Up to this Reynolds 
number Re x = 890, the important modes are the steady waves MO= 0 and the fundamental wave 
MO=l. Second observation is that beyond the secondary instability region, the disturbances 
saturate and there is no explosive growth nor any instability of the disturbances. This agrees 
with the DNS computations of Krist and Zang (1987). 

6. Conclusions 

We have developed a PSE code where the nonlinear terms are treated implicitly. The 
equations are discretized using two-point fourth order compact scheme in the y-direction and the 
backward Euler method in the x-direction. We investigated the transition process in a Blasius 
boundary layer due to fundamental type breakdown. With the 8 modes in u ; and 3 planes, we 
could compute the evolution of disturbances up to Re x =910, which is well into the strongly 
nonlinear region. The transition onset point is located around Re x =850. The comparison with 
the measurements and with the DNS computations are very good up to Re=880. 

In future, we have to repeat this computations with larger number of modes. This will require 
large memory and may not be possible to increase the modes beyond a limit. One solution is 
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to look for efficient iterative schemes to solve the coupled equations. Secondly, we have to 
perform the computations for several different cases to conclude the effectiveness of the PSE 
approach in the fully nonlinear region. 
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Figure 5. Distribution and evolution of the rms value of the axial velocity. 
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Figure 9(f). 










